This notebook analyzes the 55 Cancri system by testing the 5-planet model, with a Gaussian process (GP) to model stellar activity. After the 5-planet model is constructed and analyzed, we demonstrate our post-hoc method of estimating jitter, and compare results to our 4-planet model.
In order to run this notebook, you will need to download all of the python packages we use. The complete list of used packages can be found at the top of the "functions.py" file in the repository.
from functions import * # See functions.py for the list of used functions and packages
# Orbital parameters
f_mag = 1.0/(10.5*365.25) #This equals 1/3835.125 (Bourrier value: 1.0/3822.4)
f_rot = 1.0/(38.8) # From Bourrier et al. 2018
planet_b = 1.0/14.65314 # (Bourrier value: 1.0/14.6516)
planet_c = 1.0/44.373 #(Bourrier value: 1.0/44.3989)
planet_d = 1.0/4867 #(Bourrier value: 1.0/5574.2)
planet_e = 1.0/0.7365478 #(Bourrier value: 1.0/0.73654737)
planet_f = 1.0/260.91 #(Bourrier value: 1.0/259.88)
lunar = 1.0/29.53
planets = [planet_b, planet_c, planet_d, planet_e, planet_f]
# Read in RV and S-index
inst_Bourrier = np.loadtxt("data/Bourrier_RV.txt", dtype = str, usecols = (5), skiprows = 43)
t_Bourrier, RV_Bourrier, RVerr_Bourrier, S_Bourrier, Serr_Bourrier = np.loadtxt("data/Bourrier_RV.txt", usecols = (0, 1, 2, 3, 4), skiprows = 43, unpack = True ) # BINNED, TRENDED
t_Bourrier = t_Bourrier + 2400000.0 # Shift to JD
# Sort by time index
t_Bourrier, RV_Bourrier, RVerr_Bourrier, inst_Bourrier = sort_time(t_Bourrier, RV_Bourrier, RVerr_Bourrier, inst_Bourrier)
instruments = inst_Bourrier # Using my read-in method
fit_method = 'L-BFGS-B'
fit_options = {
'maxiter': 1000,
'maxcor': 50
}
fit_ecc = True
fap_max = 1e-3 # FAP over this value terminates the agnostic DACE search
rv_model = rv.RvModel(t_Bourrier-t_Bourrier[0], RV_Bourrier, err = term.Error(RVerr_Bourrier))
# Construct a dictionary for the RV data
rv_dict = {}
rv_dict['jd'] = t_Bourrier-t_Bourrier[0]
rv_dict['rv'] = RV_Bourrier
rv_dict['rv_err'] = RVerr_Bourrier
rv_dict['ins_name'] = instruments
# Add linear parameters
for inst in np.unique(instruments):
rv_model.add_lin(1.0*(rv_dict['ins_name']==inst), f'offset_inst_{inst}')
print("Fit linear parameters:")
rv_model.fit(method=fit_method, options=fit_options)
rv_model.show_param()
print('loglikelihood =', rv_model.loglike())
rv_err = np.sqrt(rv_model.cov.A)
plt.figure()
for inst in np.unique(instruments):
kinst = rv_dict['ins_name'] == inst
plt.errorbar(rv_dict['jd'][kinst],
rv_dict['rv'][kinst] - rv_model.get_param(f'lin.offset_inst_{inst}'),
yerr=rv_err[kinst],
fmt='.', rasterized=True, label = inst)
plt.xlabel('Time [days]')
plt.ylabel('$\Delta v$ [m/s]')
plt.legend(loc = 'lower left', fontsize = 'x-small')
plt.show()
Fit linear parameters: Parameter Value Error lin.offset_inst_HARPN 27391.0595 ± 0.0373 lin.offset_inst_HARPS 27463.0762 ± 0.0481 lin.offset_inst_HRS 28354.962 ± 0.226 lin.offset_inst_KECK -14.0313 ± 0.0434 lin.offset_inst_LICK 2.0091 ± 0.0923 lin.offset_inst_SOPHIE 27461.2180 ± 0.0982 lin.offset_inst_TULL -22557.388 ± 0.164 loglikelihood = -1555297.8373726257
RR = 1./(max(rv_dict['jd'])-rv_dict['jd'][0]) # Rayleigh Resolution, to be used in frequency analysis
df = RR/2
f_N = 1./(2*np.median(np.diff(rv_dict['jd']))) # Nyquist frequency
f = np.arange(0, f_N, df)
# This function loops over the frequencies contained within the "signals" input array, fits Keplerians to each,
# and computes the resulting residual.
def remove_kep_without_jitter(RV_DICT, RV_MODEL, signals, fit_param = True):
time, RV, rverr = RV_DICT['jd'], RV_DICT['rv'], RV_DICT['rv_err']
rv_model = RV_MODEL # Just redine for c/p reasons
N = len(time)
# Customize frequency array
RR = 1./(max(time)-time[0])
delta_f = RR/2
f_N = 1./(2*np.median(np.diff(rv_dict['jd']))) # Nyquist frequency
f = np.arange(0, f_N, delta_f)
N_signals = len(signals)
loglike = np.zeros(N_signals)
LR = np.zeros(N_signals - 1)
BIC = np.zeros(N_signals - 1)
p_val = np.zeros(N_signals - 1)
df = 5 # Degrees of freedom in Keplerian model
alpha = 0.01 # Critical p-value
LR_cutoff = chi2.ppf(1-alpha, df) # Cutoff for LR hypothesis testing
# Loop through, removing planets
for passthrough in range(N_signals):
res = rv_model.residuals() # Update residual
rv_err = np.sqrt(rv_model.cov.A)
# Plot the RV residual scatterplot
plot_residual(RV_DICT, res, rv_model, passthrough)
fls, pls, GLS_window, percentiles = GLS_calc(time, res, bootstrap = True)
# Plot the GLS periodogram of the residual
plot_GLS(fls, pls, percentiles, passthrough)
print('Removing orbit of P={} d \n'.format(1.0/signals[passthrough]))
rv_model.add_keplerian_from_period((1.0/signals[passthrough]), fit = True)
rv_model.set_keplerian_param(f'{rv_model.nkep-1}', param=['P', 'la0', 'K', 'e', 'w'])
if not fit_ecc:
rv_model.set_param(np.zeros(2), rv_model.fit_param[-2:])
rv_model.fit_param = rv_model.fit_param[:-2]
# Global fit of the model
if fit_param:
print("Fitting model parameters: ")
rv_model.fit(method=fit_method, options=fit_options)
rv_model.show_param()
print('loglikelihood = {}\n'.format(rv_model.loglike()))
#print("Residual RMS: ", np.sqrt(np.mean(res**2)))
#print("Residual Chi^2: ", rv_model.chi2())
loglike[passthrough] = rv_model.loglike()
# After the final iteration of the loop, plot the final resiudal scatter plot and GLS
res = rv_model.residuals() # Final Update of residual
rv_err = np.sqrt(rv_model.cov.A)
fls, pls, GLS_window, percentiles = GLS_calc(time, res, bootstrap = True)
plot_residual(RV_DICT, res, rv_model, (passthrough+1))
plot_GLS(fls, pls, percentiles, passthrough+1)
# Compute Likelihood Ratio test for the nested planetary models
for i in range(len(LR)):
LR[i] = 2*(loglike[i+1] - loglike[i]) # Order is less restrictive model minus more restrictive (more params - less params)
p_val[i] = chi2.sf(LR[i], df) # 1 - cdf
BIC[i] = len(rv_model.get_param())*np.log(len(RV)) - 2*loglike[i]
print("LR cutoff for significance: ", LR_cutoff)
print("\t ~~ 55Cnc No Jitter ~~", "\nloglike : ", loglike, "\nLR = ", LR, "\n")
print("p-values: ", p_val, "\n")
print("Bayesian information criterion: ", BIC)
return rv_model, res, loglike, LR, p_val
Let's compute the 5-Keplerian model. 5 instances of bootstrapped FALs in the periodograms will take some time. $P = 391 = N/2$ days is a spectral window artifact.
model_5pGP, res, loglike, LR, p_val = remove_kep_without_jitter(rv_dict, rv_model, [planet_b, planet_c, planet_f, planet_e, planet_d], fit_param = True)
0.0 0.1 0.2 0.4 0.8 Bootstrap took 100.73886609077454 s for 10000 iterations
Removing orbit of P=14.653140000000002 d Fitting model parameters: Parameter Value Error lin.offset_inst_HARPN 27440.5723 ± 0.0503 lin.offset_inst_HARPS 27479.3079 ± 0.0608 lin.offset_inst_HRS 28357.591 ± 0.226 lin.offset_inst_KECK -23.7175 ± 0.0448 lin.offset_inst_LICK 4.7905 ± 0.0926 lin.offset_inst_SOPHIE 27449.1494 ± 0.0995 lin.offset_inst_TULL -22558.505 ± 0.164 kep.0.P 14.6523452 ± 0.0000153 kep.0.la0 [deg] 336.564 ± 0.204 kep.0.K 69.4968 ± 0.0449 kep.0.e 0.023049 ± 0.000541 kep.0.w [deg] 203.78 ± 1.33 loglikelihood = -262917.29213645123
0.0 0.1 0.2 0.4 0.8 Bootstrap took 98.70048522949219 s for 10000 iterations
Removing orbit of P=44.373 d Fitting model parameters: Parameter Value Error lin.offset_inst_HARPN 27444.636 ± 0.119 lin.offset_inst_HARPS 27478.4958 ± 0.0821 lin.offset_inst_HRS 28356.379 ± 0.231 lin.offset_inst_KECK -23.2234 ± 0.0575 lin.offset_inst_LICK 4.520 ± 0.101 lin.offset_inst_SOPHIE 27446.327 ± 0.101 lin.offset_inst_TULL -22557.963 ± 0.164 kep.0.P 14.6517643 ± 0.0000199 kep.0.la0 [deg] 329.838 ± 0.259 kep.0.K 68.7014 ± 0.0481 kep.0.e 0.018235 ± 0.000569 kep.0.w [deg] 114.89 ± 2.37 kep.1.P 44.36852 ± 0.00368 kep.1.la0 [deg] 318.35 ± 5.45 kep.1.K 10.7063 ± 0.0665 kep.1.e 0.3585 ± 0.0138 kep.1.w [deg] 124.503 ± 0.781 loglikelihood = -221290.61866871596
0.0 0.1 0.2 0.4 0.8 Bootstrap took 96.31640338897705 s for 10000 iterations
Removing orbit of P=260.91 d Fitting model parameters:
/home/justin/anaconda3/envs/freshstart/lib/python3.10/site-packages/kepmodel/timeseries.py:716: RuntimeWarning: invalid value encountered in sqrt error = np.sqrt(-np.diag(np.linalg.inv(hess)))
Parameter Value Error lin.offset_inst_HARPN 27447.2918 ± 0.0632 lin.offset_inst_HARPS 27475.3913 ± 0.0694 lin.offset_inst_HRS 28357.348 ± 0.227 lin.offset_inst_KECK -22.5203 ± 0.0437 lin.offset_inst_LICK 4.0621 ± 0.0930 lin.offset_inst_SOPHIE 27443.449 ± 0.102 lin.offset_inst_TULL -22557.678 ± 0.164 kep.0.P 14.6516854 ± 0.0000158 kep.0.la0 [deg] 327.377 ± 0.207 kep.0.K 68.8336 ± 0.0462 kep.0.e 0.024540 ± 0.000552 kep.0.w [deg] 128.31 ± 1.33 kep.1.P 44.46008 ± 0.00104 kep.1.la0 [deg] 94.44 ± 1.52 kep.1.K 10.2263 ± 0.0404 kep.1.e 0.18287 ± 0.00365 kep.1.w [deg] 141.14 ± 1.28 kep.2.P 261.29028 ± 0.00299 kep.2.la0 [deg] 343.040847 ± nan kep.2.K 25.219214 ± nan kep.2.e 0.950000 ± nan kep.2.w [deg] 123.326085 ± nan loglikelihood = -205447.7204552926
0.0 0.1 0.2 0.4 0.8 Bootstrap took 99.21420550346375 s for 10000 iterations
Removing orbit of P=0.7365478 d Fitting model parameters: Parameter Value Error lin.offset_inst_HARPN 27445.4842 ± 0.0694 lin.offset_inst_HARPS 27473.4502 ± 0.0766 lin.offset_inst_HRS 28357.368 ± 0.227 lin.offset_inst_KECK -22.4317 ± 0.0439 lin.offset_inst_LICK 3.7851 ± 0.0933 lin.offset_inst_SOPHIE 27441.565 ± 0.107 lin.offset_inst_TULL -22557.928 ± 0.165 kep.0.P 14.6517655 ± 0.0000159 kep.0.la0 [deg] 328.531 ± 0.209 kep.0.K 69.0084 ± 0.0464 kep.0.e 0.025377 ± 0.000551 kep.0.w [deg] 129.28 ± 1.28 kep.1.P 44.45957 ± 0.00103 kep.1.la0 [deg] 93.70 ± 1.50 kep.1.K 10.3062 ± 0.0403 kep.1.e 0.17643 ± 0.00360 kep.1.w [deg] 140.88 ± 1.30 kep.2.P 261.29427 ± 0.00238 kep.2.la0 [deg] 346.010585 ± nan kep.2.K 25.862268 ± nan kep.2.e 0.950000 ± nan kep.2.w [deg] 126.177080 ± nan kep.3.P 0.736549 ± nan kep.3.la0 [deg] 124.364 ± 0.376 kep.3.K 26.564 ± 0.927 kep.3.e 0.95000 ± 0.00223 kep.3.w [deg] 167.287 ± 0.383 loglikelihood = -201799.83301960386
0.0 0.1 0.2 0.4 0.8 Bootstrap took 95.99604082107544 s for 10000 iterations
Removing orbit of P=4867.0 d Fitting model parameters: Parameter Value Error lin.offset_inst_HARPN 27459.218 ± 0.158 lin.offset_inst_HARPS 27460.000 ± 0.150 lin.offset_inst_HRS 28395.260 ± 0.259 lin.offset_inst_KECK -42.6174 ± 0.0762 lin.offset_inst_LICK 0.7857 ± 0.0959 lin.offset_inst_SOPHIE 27430.713 ± 0.164 lin.offset_inst_TULL -22576.028 ± 0.181 kep.0.P 14.6513954 ± 0.0000153 kep.0.la0 [deg] 324.729 ± 0.203 kep.0.K 71.0776 ± 0.0470 kep.0.e 0.006467 ± 0.000549 kep.0.w [deg] 90.18 ± 5.16 kep.1.P 44.42064 ± 0.00103 kep.1.la0 [deg] 37.29 ± 1.48 kep.1.K 10.1576 ± 0.0377 kep.1.e 0.04614 ± 0.00404 kep.1.w [deg] 75.66 ± 4.91 kep.2.P 261.0661 ± 0.0669 kep.2.la0 [deg] 335.06 ± 2.78 kep.2.K 4.7572 ± 0.0409 kep.2.e 0.000000000000 ± 0.000000000278 kep.2.w [deg] 116.471020 ± nan kep.3.P 0.736552 ± nan kep.3.la0 [deg] 128.651563 ± nan kep.3.K 6.2004 ± 0.0443 kep.3.e 0.05032 ± 0.00777 kep.3.w [deg] 139.33 ± 4.94 kep.4.P 5187.12 ± 6.10 kep.4.la0 [deg] 200.478 ± 0.563 kep.4.K 43.0793 ± 0.0805 kep.4.e 0.11611 ± 0.00141 kep.4.w [deg] 287.305 ± 0.894 loglikelihood = -12207.151009076284 0.0 0.1 0.2 0.4 0.8 Bootstrap took 95.92229557037354 s for 10000 iterations
LR cutoff for significance: 15.08627246938899 ~~ 55Cnc No Jitter ~~ loglike : [-262917.29213645 -221290.61866872 -205447.72045529 -201799.8330196 -12207.15100908] LR = [ 83253.34693547 31685.79642685 7295.77487138 379185.36402106] p-values: [0. 0. 0. 0.] Bayesian information criterion: [526047.7636246 442794.41668913 411108.62026228 403812.84539091]
time = rv_dict['jd'] - rv_dict['jd'][0]; RV = RV_Bourrier; RVerr = RVerr_Bourrier
unif_time = np.linspace(time[0], np.max(time), round(len(time)/3))
kep_model = model_5pGP
N_param = len(rv_model.get_param())
metric_guess = 17.0
gamma_guess = 1.0
period_guess = np.exp(8.55)
# Set reasonable boundaries for each hyperparameter
# Parameter order: ln(amplitude), ln(metric^2), gamma, ln(period)
GP_lower_bounds = [0.5, np.log(2000.0**2), 0.01, np.log(1500.0)]
GP_upper_bounds = [5.8, 22.0, 10.0, 8.8]
model_values = np.array(kep_model.get_param())
model_errors = np.array(kep_model.get_param_error()[1])
# Fix the nan in uncertanties to Bourrier et al (2018) value
model_errors[N_param-10] = 1.3*10**(-6)
# Vastly increasing parameter error bars to allow GP to explore more space
# (This only affects instrumental offset bounds; kep param bounds fixed below)
kep_lower_bounds = model_values - 1000.0*model_errors
kep_upper_bounds = model_values + 1000.0*model_errors
for p in range(5):
kep_lower_bounds[7 + int(5*p) + 1] = 0.0 # Mean longitude
kep_lower_bounds[7 + int(5*p) + 3] = 0.0 # eccentricity
kep_lower_bounds[7 + int(5*p) + 4] = 0.0 # argument of periastron
kep_upper_bounds[7 + int(5*p) + 1] = 360.0 # Mean longitude
kep_upper_bounds[7 + int(5*p) + 3] = 0.7 # eccentricity
kep_upper_bounds[7 + int(5*p) + 4] = 360.0 # argument of periastron
# Fix amplitude bounds
kep_lower_bounds[16-7] = 60.0; kep_upper_bounds[16-7] = 80.0 # Amplitude of planet b
kep_lower_bounds[21-7] = 2.0; kep_upper_bounds[21-7] = 20.0 # Amplitude of planet c
kep_lower_bounds[26-7] = 2.0; kep_upper_bounds[26-7] = 20.0 # Amplitude of planet f
kep_lower_bounds[26-9] = 200.0; kep_upper_bounds[26-9] = 300.0 # Period of planet f
kep_lower_bounds[31-7] = 2.0; kep_upper_bounds[31-7] = 20.0 # Amplitude of planet e
kep_lower_bounds[N_param - 3] = 10.0; kep_upper_bounds[N_param - 3] = 50.0 # Amplitude of planet d
kep_lower_bounds[N_param - 5] = 1500.0; kep_upper_bounds[N_param - 5] = 6500.0 # Period of planet d
par_bounds = Bounds(combine_kep_GP_params(kep_lower_bounds, GP_lower_bounds), combine_kep_GP_params(kep_upper_bounds, GP_upper_bounds))
k_exp2 = np.std(kep_model.residuals()) * kernels.ExpSquaredKernel(metric = metric_guess)
k_per = kernels.ExpSine2Kernel(gamma = gamma_guess, log_period = np.log(period_guess))
k_mag = k_exp2 * k_per
# The below guess parameters are prior Nelder-Mead results
guess_pars_no_jitter = [ 2.74550055e+04, 2.74691022e+04, 2.83974958e+04, -3.94343675e+01,
4.24743467e+00, 2.74390579e+04, -2.25703040e+04, 1.46515548e+01,
5.70595877e+00, 7.12624543e+01, 7.69216943e-04, 3.91169320e+00,
4.44030098e+01, 1.79906272e-01, 9.72853639e+00, 3.43073964e-02,
6.52905190e+00, 2.59874924e+02, 5.12012594e+00, 5.22035182e+00,
2.36240133e-01, 5.34052387e+00, 7.36548081e-01, 1.94569550e+00,
6.01660310e+00, 3.27734161e-02, 2.07926296e+00, 5.62084890e+03,
3.58342626e+00, 2.45989628e+01, 3.82921437e-01, 4.92326155e+00,
3.47449588e+00, 1.64768058e+01, 3.03392484e+00, 8.42431518e+00]
gp_irregular = george.GP(k_mag, mean = np.mean(kep_model.residuals()), fit_kernel = True)
gp_irregular.set_parameter_vector(guess_pars_no_jitter[len(guess_pars_no_jitter)-4:])
gp_irregular.compute(time, RVerr)
Below is my previous best fit GP for the 5 planet model
for i in range(len(kep_lower_bounds)):
print("Min | Initial Value | Max : ", kep_lower_bounds[i], "|", guess_pars_no_jitter[i], "|", kep_upper_bounds[i])
Min | Initial Value | Max : 27301.09865381575 | 27455.0055 | 27617.33770768819 Min | Initial Value | Max : 27310.0282411373 | 27469.1022 | 27609.971162289065 Min | Initial Value | Max : 28136.17475011404 | 28397.4958 | 28654.34443994912 Min | Initial Value | Max : -118.78867099043451 | -39.4343675 | 33.55377966284321 Min | Initial Value | Max : -95.11071606518351 | 4.24743467 | 96.68219007535328 Min | Initial Value | Max : 27266.79707516103 | 27439.0579 | 27594.62861882425 Min | Initial Value | Max : -22757.471279995272 | -22570.304 | -22394.583728498266 Min | Initial Value | Max : 14.636069880612531 | 14.6515548 | 14.66672087347186 Min | Initial Value | Max : 0.0 | 5.70595877 | 360.0 Min | Initial Value | Max : 60.0 | 71.2624543 | 80.0 Min | Initial Value | Max : 0.0 | 0.000769216943 | 0.7 Min | Initial Value | Max : 0.0 | 3.9116932 | 360.0 Min | Initial Value | Max : 43.3953213655894 | 44.4030098 | 45.44595766323891 Min | Initial Value | Max : 0.0 | 0.179906272 | 360.0 Min | Initial Value | Max : 2.0 | 9.72853639 | 20.0 Min | Initial Value | Max : 0.0 | 0.0343073964 | 0.7 Min | Initial Value | Max : 0.0 | 6.5290519 | 360.0 Min | Initial Value | Max : 200.0 | 259.874924 | 300.0 Min | Initial Value | Max : 0.0 | 5.12012594 | 360.0 Min | Initial Value | Max : 2.0 | 5.22035182 | 20.0 Min | Initial Value | Max : 0.0 | 0.236240133 | 0.7 Min | Initial Value | Max : 0.0 | 5.34052387 | 360.0 Min | Initial Value | Max : 0.7352520414971537 | 0.736548081 | 0.7378520414971537 Min | Initial Value | Max : 0.0 | 1.9456955 | 360.0 Min | Initial Value | Max : 2.0 | 6.0166031 | 20.0 Min | Initial Value | Max : 0.0 | 0.0327734161 | 0.7 Min | Initial Value | Max : 0.0 | 2.07926296 | 360.0 Min | Initial Value | Max : 1500.0 | 5620.8489 | 6500.0 Min | Initial Value | Max : 0.0 | 3.58342626 | 360.0 Min | Initial Value | Max : 10.0 | 24.5989628 | 50.0 Min | Initial Value | Max : 0.0 | 0.382921437 | 0.7 Min | Initial Value | Max : 0.0 | 4.92326155 | 360.0
kep_model.set_param(guess_pars_no_jitter[0:len(guess_pars_no_jitter)-4])
gp_irregular.set_parameter_vector(guess_pars_no_jitter[len(guess_pars_no_jitter)-4:])
gp_irregular.compute(time, RVerr)
t_GP_TEST, GP_fit_TEST = plot_GP(gp_irregular, gp_irregular.get_parameter_vector(), time, kep_model.residuals(), RVerr)
To model stellar activity, we construct a quasiperiodic GP kernel as informed by Equation 27 from Rajpaul et al. (2015): \begin{equation} k_{\text{QP}}(t_i, t_j) = K^2 \text{exp}\Bigg\{ - \frac{ \sin^2{[\pi(t_i - t_j)/P]}}{2 \lambda_p^2} - \frac{(t_i - t_j)^2}{2 \lambda_e^2} \Bigg\} , \label{eq:GP_ker} \end{equation} where $K$ is the amplitude, $P$ is the period of oscillation, $\lambda_p$ is the dimensionless length scale of periodic variation (often called a correlation time scale), and $\lambda_e$ is an evolutionary time scale in the units of $t$. The first term captures periodic variation through the non-stationary exponential-sine-squared kernel and amplitude modulation though the stationary exponential-squared kernel.
Our model fitting is through a curvature-penalized likelihood function $\mathcal{L}(\boldsymbol{\theta})$ given input parameters $\boldsymbol{\theta}$: \begin{equation} \mathcal{L}(\boldsymbol{\theta}) = -\frac{1}{2} \sum^{N-1}_{n = 0} \frac{(\Delta y_n - \tau_n)^2}{\sigma_{y, n}^2} - \frac{1}{2} \alpha \sum^{M}_{m = 1} \Bigg( \frac{\tau_{m+1} - 2\tau_m + \tau_{m-1}}{\Delta t_{\text{GP}}^2} \Bigg)^2 \end{equation} where $\Delta y_n$ is the Keplerian residual and $\tau_n$ is the stellar activity induced Doppler shift predicted by $k_{\text{QP}}$ at time $t_n$. $\alpha$ dictates the degree to which GP curvature is penalized.
You can play around with Nelder-Mead optimization using the code below. We ultimately ran both Nelder-Mead and MCMC code on the University of Delaware Caviness computing cluster. The MCMC script can be found in the main repository branch, and is appropriately labeled.
par_bounds = Bounds(combine_kep_GP_params(kep_lower_bounds, GP_lower_bounds), combine_kep_GP_params(kep_upper_bounds, GP_upper_bounds))
alpha = 3.0e9
def penalized_NLL(all_params):
kep_params, GP_params = split_kep_gp(all_params)
kep_model.set_param(kep_params)
kepres = kep_model.residuals() # Update residual
gp_irregular.set_parameter_vector(GP_params)
# Compute mean GP prediction
GP_pred_unif, pred_var_unif = gp_irregular.predict(kepres, unif_time, return_var = True)
GP_pred_irregular, pred_var_irregular = gp_irregular.predict(kepres, time, return_var = True)
deriv_sum = (second_deriv(unif_time, GP_pred_unif)**2).sum()
obj = -0.5*((kepres-GP_pred_irregular)**2/RVerr_Bourrier**2).sum() - 0.5*alpha*deriv_sum
return -obj if np.isfinite(obj) else 1e25
def penalized_NLL_PRINTMODE(all_params): # Same as above, just prints out the values of the two terms in the obj func
kep_params, GP_params = split_kep_gp(all_params)
kep_model.set_param(kep_params)
kepres = kep_model.residuals() # Update residual
gp_irregular.set_parameter_vector(GP_params)
# Compute mean GP prediction
GP_pred_unif, pred_var_unif = gp_irregular.predict(kepres, unif_time, return_var = True)
GP_pred_irregular, pred_var_irregular = gp_irregular.predict(kepres, time, return_var = True)
deriv_sum = (second_deriv(unif_time, GP_pred_unif)**2).sum()
obj = -0.5*((kepres-GP_pred_irregular)**2/RVerr_Bourrier**2).sum() - 0.5*alpha*deriv_sum
print("LSQ term = ", -0.5*((kepres-GP_pred_irregular)**2/RVerr_Bourrier**2).sum())
print("Curvature term = ", - 0.5*alpha*deriv_sum)
print("alpha = ", alpha)
return -obj if np.isfinite(obj) else 1e25
print("\nBeginning O(p) = {0:.2f}".format(penalized_NLL(guess_pars_no_jitter)))
print("LAST RUN GP params: ")
print("Fitted GP period [days]: ", np.exp(guess_pars_no_jitter[len(guess_pars_no_jitter)-1]))
print("GP decor. timescale [days]: ", np.sqrt(np.exp(guess_pars_no_jitter[len(guess_pars_no_jitter)-3])) )
print("Fitted GP amplitude [m/s]: ", np.exp(guess_pars_no_jitter[len(guess_pars_no_jitter)-4]))
optimize = False
if optimize:
print(penalized_NLL_PRINTMODE(guess_pars_no_jitter))
start = timer.time()
result = minimize(penalized_NLL, guess_pars_no_jitter, method="Nelder-Mead", bounds = par_bounds,
options = {'disp': True, 'maxfev': 2000})
end = timer.time()
print("Optimization took ", (end-start)/60.0, "min")
print("\nFinal O(p) = {0:.2f}".format(penalized_NLL(result.x)))
print(penalized_NLL_PRINTMODE(result.x))
print("=========Optimized Parameters=========\n", result.x)
print("Fitted GP period [days]: ", np.exp(result.x[len(result.x)-1]))
print("GP decor. timescale [days]: ", np.sqrt(np.exp(result.x[len(result.x)-3])) )
print("Fitted GP amplitude [m/s]: ", np.exp(result.x[len(result.x)-4]))
print("Success?: ", result.success)
#print(result.x)
Beginning O(p) = 7094.93 LAST RUN GP params: Fitted GP period [days]: 4556.523311404423 GP decor. timescale [days]: 3783.4928638165625 Fitted GP amplitude [m/s]: 32.28155065863546
We computed MCMC results using the University of Delaware Caviness cluster, and simply read in the resulting chains for further analysis. See "55Cnc_5pGP_MCMC_alpha_3e9_final.py" for the MCMC code.
chains_3e9_newamp = np.loadtxt('caviness/55Cnc_5pGP_MCMC_chains_alpha_3000000000_newamp.txt')
Let's investigate the MCMC results for the curvature penalty we end up using for this model: $\alpha = 3\cdot10^9$.
nsamples = 400000
alpha = 3.0e9
samples = chains_3e9_newamp.reshape((nsamples+1, 5*5+7+4)) # 25 kep params, 7 inst offsets, 4 GP params
burn = nsamples//3
R = acf.acf(samples[burn:, :])
tau = np.arange(R.shape[0])
plt.figure()
plt.plot(tau[1:], R[1:])
plt.xscale('log')
plt.xlim(1, tau[-1])
plt.xlabel('Lag')
plt.ylabel('ACF')
iat = acf.iat(R=R)
iat_max = np.max(acf.iat(R=R))
ESS = R.shape[0]/iat
print('Maximum IAT:', iat_max) # Time needed for chain to forget where it started
print('Effective Number of Samples per param:', ESS) # N/iat
print("Mean ESS: ", np.mean(ESS))
print("Minimum ESS: ", np.min(ESS))
Maximum IAT: 551.7498454833271 Effective Number of Samples per param: [ 767.20123253 764.07366752 765.93773199 764.36961417 766.11266101 765.33234037 765.26614052 1046.63458387 1040.81717197 1100.8658704 931.41561699 483.31323005 1161.1975855 1191.91847251 1045.41438541 1191.88294271 1094.02117627 1075.31828504 1095.86914604 1068.34803625 1059.73412096 925.08351641 1123.52589215 1093.27117724 1067.42960853 1031.01223004 1074.23034107 1027.16139579 1032.37848125 960.58884413 1118.09163812 1081.04230016 834.83085475 735.73379505 888.82336795 907.07192715] Mean ESS: 967.9255383856726 Minimum ESS: 483.3132300497549
# Set up labels for corner plot
corner_labels = ["HARPN "+"$\Delta v$", "HARPS "+"$\Delta v$", "HRS "+"$\Delta v$", "KECK "+"$\Delta v$", "LICK "+"$\Delta v$", "SOPHIE "+"$\Delta v$", "TULL "+"$\Delta v$"]
for p in range(5):
if p==0:
corner_labels.append("(B) "+"$P$")
corner_labels.append("(B) "+"$L_0$")
corner_labels.append("(B) "+"$K$")
corner_labels.append("(B) "+"$e$")
corner_labels.append("(B) "+"$\omega$")
elif p==1:
corner_labels.append("(C) "+"$P$")
corner_labels.append("(C) "+"$L_0$")
corner_labels.append("(C) "+"$K$")
corner_labels.append("(C) "+"$e$")
corner_labels.append("(C) "+"$\omega$")
elif p==2:
corner_labels.append("(F) "+"$P$")
corner_labels.append("(F) "+"$L_0$")
corner_labels.append("(F) "+"$K$")
corner_labels.append("(F) "+"$e$")
corner_labels.append("(F) "+"$\omega$")
elif p==3:
corner_labels.append("(E) "+"$P$")
corner_labels.append("(E) "+"$L_0$")
corner_labels.append("(E) "+"$K$")
corner_labels.append("(E) "+"$e$")
corner_labels.append("(E) "+"$\omega$")
elif p==4:
corner_labels.append("(D) "+"$P$")
corner_labels.append("(D) "+"$L_0$")
corner_labels.append("(D) "+"$K$")
corner_labels.append("(D) "+"$e$")
corner_labels.append("(D) "+"$\omega$")
corner_labels.append("$\ln{(K)}$")
corner_labels.append(" $\ln{(\lambda_{e}^2)}$")
corner_labels.append("$1/2\lambda_{p}^2$")
corner_labels.append(" $P$")
# Corner plot of only the instrumentation effects + a couple planets
print("Burn-in: ", burn)
corner.corner(
samples[burn:, 0:17], #samples[nsamples//4:, 19:],
labels=kep_model.fit_param[0:17], #corner_labels[0:17], #kep_model.fit_param[0:17],
quantiles=[0.159, 0.5, 0.841],
show_titles=True)
plt.show()
plt.close()
Burn-in: 133333
label_dict = {}; label_dict["fontsize"] = "xx-large"
title_kwargs = {}; title_kwargs["fontsize"] = "xx-large"
corner.corner(
samples[burn:, 17:],
labels=corner_labels[17:],
quantiles=[0.159, 0.5, 0.841],
show_titles=True, max_n_ticks=3, label_kwargs = label_dict, title_kwargs = title_kwargs)
plt.show()
plt.close()